Aluminum (Al3+) toxicity is a typical abiotic stress that severely limits crop production in acidic soils. In this study, an RIL (recombinant inbred line, F12) population derived from the cross of Zhonghuang 24 (ZH 24) and Huaxia 3 (HX 3) (160 lines) was tested using hydroponic cultivation. Relative root elongation (RRE) and apical Al3+ content (AAC) were evaluated for each line, and a significant negative correlation was detected between the two indicators. Based on a high-density genetic linkage map, the phenotypic data were used to identify quantitative trait loci (QTLs) associated with these traits. With composite interval mapping (CIM) of the linkage map, five QTLs that explained 39.65% of RRE and AAC variation were detected on chromosomes (Chrs) Gm04, Gm16, Gm17 and Gm19. Two new QTLs, qRRE_04 and qAAC_04, were located on the same region of bin93-bin94 on Chr Gm04, which explained 7.09% and 8.98% phenotypic variation, respectively. Furthermore, the results of the expression analysis of candidate genes in the five genetic regions of the QTLs showed that six genes (Glyma.04g218700, Glyma.04g212800, Glyma.04g213300, Glyma.04g217400, Glyma.04g216100 and Glyma.04g220600) exhibited significant differential expression between the Al3+ treatment and the control of two parents. The results of qRT-PCR analysis indicated that Glyma.04g218700 was upregulated by Al3+ treatment with the hundreds-fold increased expression level and may be a candidate gene with potential roles in the response to aluminum stress. Therefore, our efforts will enable future functional analysis of candidate genes and will contribute to the strategies for improvement of aluminum tolerance in soybean.

Citation: Wang X, Cheng Y, Yang C, Yang C, Mu Y, Xia Q, et al. (2019) QTL mapping for aluminum tolerance in RIL population of soybean (Glycine max L.) by RAD sequencing. PLoS ONE 14(10): e0223674. https://doi.org/10.1371/journal.pone.0223674

Editor: Harsh Raman, New South Wales Department of Primary Industries, AUSTRALIA

Received: May 12, 2019; Accepted: September 25, 2019; Published: October 29, 2019

Copyright: © 2019 Wang et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: All relevant data are within the paper and its Supporting Information files.

Funding: This work was supported by The Project of Science and Technology of Guangzhou (201804020015)-(Hai Nian, prepared the soybean materials used in this study); The Major Project of New Varieties Cultivation of Genetically Modified Organisms (2016ZX08004002-007)-(Qibin Ma, designed the research program; The Project of Molecular Design Breeding for Major Economic Crops (2016YFD0101901)-(Cunyi Yang, designed the research program; The Project of the National Natural Sciences Foundation of China (31771816)-(Qibin Ma); The Project of the National Key R and D Program of China (2017YFD0101500)-(Hai Nian); The Project of the China Agricultural Research System (CARS-04-PS09)-(Hai Nian.

Competing interests: The authors have declared that no competing interests exist.

Abbreviations: RIL, Recombinant inbred line; QTL, Quantitative trait locus; CIM, Composite interval mapping; RRE, Relative root elongation; AAC, Apical Al3+ content; Chr, Chromosome; LOD, Log-likelihood; RFLP, Restriction fragment length polymorphism; AFLP, Amplified fragment length polymorphism; SSR, Simple sequence repeats; SNP, Single-nucleotide polymorphism; RAD-seq, Restriction-site-associated DNA sequencing; NGS, Next-generation sequencing

Aluminum (Al3+) toxicity is one of the major factors affecting crop production on acidic soils worldwide [1, 2]. When the soil pH decreases to values less than 5.0, Al is solubilized as the phytotoxic Al3+, which has a pernicious effect on crops. It was found that root elongation can be inhibited in seconds at micromolar concentrations of Al3+ [3]. The primary location of Al3+ toxicity is at the root tip where Al3+ binds to the cell wall [4]. Changes in some components of the cell wall lead to a limited capacity of damaged roots for absorption of sufficient water and nutrients from soil [2]. Additionally, the damaged root impeded the growth of shoot, and eventually reduced the yield of crops. Soybean is one of the most important crops in the subtropical zone and is also damaged by Al3+ toxicity in acidic soil. Hence, investigation of the traits associated with Al3+ toxicity via a combination of identified soybean germplasms and sequencing technology is of great significance.

It is well known that two types of mechanisms of Al3+ resistance in soybean are involved in the exclusion of Al3+ from the root apex (external exclusion) or in conferring tolerance to Al3+ in the plant symplast (internal tolerance) [5]. The mechanism of external exclusion involves excretion of organic acids to chelate Al3+ from the root cells, increasing the rhizospheric pH and external exclusion of border cells [6, 7]. However, the mechanism of internal tolerance depends on chelation of organic acids and segregation of Al3+ in vacuoles. Antioxidant metabolism as well as hormone signal transduction also contribute to aluminum tolerance [8, 9].

Aluminum tolerance of soybean is a complex quantitative trait with substantial genetic variation [10]. Studies on the genetic architecture of soybean aluminum tolerance remain challenging due to the interactions of environments and genotypes. Conventional breeding has relied on the selection of highly Al3+-tolerant cultivars for crop improvement, but this method is costly and time consuming [11]. In recent years, genome-wide association study (GWAS) and QTL mapping are commonly used to map genetic markers associated with quantitative traits. GWAS analysis generally involves the natural populations to detect the correlation between genetic polymorphism and phenotypic variation by statistical methods based on the linkage disequilibrium. A number of important GWAS loci and candidate genes for Al3+-tolerant traits have been reported over the recent decade [12]. Meanwhile, the strategy of QTL mapping has provided an improved understanding of the genetic architectures of complex traits, which has accelerated crop improvement [13]. Accordingly, extensive efforts have been directed at QTL mapping for aluminum tolerance in Arabidopsis thaliana [14] and several crops, including rice [15, 16], wheat [17, 18], barley [19], maize [20], soybean [21] and alfalfa [22].

In soybean, some QTLs of aluminum tolerance have been identified using populations from different genetic backgrounds, for which the traits of root elongation were often used to represent aluminum tolerance. In the early 2000s, a genetic linkage map containing 155 restricted fragment length polymorphism (RFLP) markers was constructed using the population derived from Young × PI 416937. Bianchihall et al. [23] detected the genetic basis of Al-tolerant traits in soybean using the map and indicated five independent RFLP markers associated with root elongation. Qi et al. [20] and Korir et al. [24] focused on the progenies of Kefeng No.1 × Nannong 1138–2 and used genetic linkage map with RFLP and simple sequence repeats (SSR) markers to detect one major and two minor QTLs for aluminum tolerance. In general, the explorations of QTL mapping indicated that approximately two to five dominant loci controlled the variation in Al-tolerance levels [25, 26].

However, the traditional molecular markers, including RFLP, SSR and amplified fragment length polymorphism (AFLP), exhibited low density and uneven distribution throughout the genome [27]. QTL mapping of complex quantitative traits such as aluminum tolerance on soybean remains elusive due to the limited efficiency and accuracy of QTL positioning. In recent years, single-nucleotide polymorphism (SNP) markers have emerged with the assistance of high-throughput sequencing technology and have been mapped across plant genomes with high density and relatively even distributions, thereby improving the accuracy of QTL mapping. Over the last few years, high-density genetic maps have been constructed using recombination bins as markers [28]. Restriction-site-associated DNA sequencing (RAD-seq) [29], one of the next-generation sequencing (NGS) methods [30], has been effectively used for high-density SNP marker discovery and QTL analysis [31, 32]. In barley and wheat, high-density genetic maps have been established using RAD-seq technology with hundreds of thousands of SNP markers as well as other polymorphic markers [33]. Abdel-Haleem et al. improved the linkage map using the progenies derived from the cross of Young and PI416937 and further developed Glyma08g42400-SNP as a major QTL to be used for marker-assisted selection of aluminum tolerance [34]. Recently, a high-density genetic linkage map based on RAD-seq technology was constructed to map QTLs for both yield-related and quality traits [35, 36]. The genetic maps with ultrahigh density for the complex polyploid crops with DNA markers indicate that RAD-seq technology can be practically applied to identify the genetic basis of complex quantitative traits.

The objectives of the present study were to develop a high-density genetic map using bin markers with RAD-seq technology to identify QTLs for the traits of aluminum tolerance in the F12 RIL population derived from the cross of Zhonghuang 24 (ZH 24) and Huaxia 3 (HX 3) and to analyse candidate genes that may influence aluminum tolerance using Gene Ontology (GO) enrichment analysis.

An RIL population with 160 lines of the F12 generation derived from a cross between ZH 24 (female parent) and HX 3 (male parent) was used in the current study. ZH 24 is an Al3+-sensitive cultivar derived from Fendou 31 × Zhongdou 19, while HX 3 is an Al3+-tolerant cultivar derived from Guizao 1 × BRSMG68 (a high-yield Brazilian cultivar) [35]. All the F12 lines of the RIL population and their parents were provided by the Guangdong Subcenter of National Center for Soybean Improvement, South China Agricultural University.

A preliminary test was designed to determine the appropriate concentration of Al3+ and Al3+ treatment for hydroponic cultivation. The two parents and five randomly selected lines (L10, L70, L154, L206, and L245) were used to identify Al3+ tolerance with RRE as a detection index. The concentrations of AlCl3 (0.5 mM CaCl2, pH 4.5) were set as 0, 5, 15, 20, 25 and 30 μΜ. The RRE of each line and parents was measured by imaging analysis during successive treatment periods of 24 h, 48 h and 72 h. The Al3+ concentration and duration that provided the widest separation among these lines were chosen for screening the RIL population.

The phenotype of the RIL population was estimated by the RRE and AAC after hydroponic cultivation along with the parents. For each line as well as the two parents, the hydroponic experiments were carried out with three replications. For each replication, 6 seedlings with nearly the same root length (approximately 8 cm) were fixed using sponge in the holes of foam floating plate in plastic containers either with or without AlCl3 treatment (0.5 mM CaCl2, pH 4.5). The average values of phenotypic data for RRE and AAC were used for mapping and identifying QTLs for aluminum tolerance.

A total of 80–100 plump seeds of each line and the parents were germinated in sterilized vermiculite for three days at 26°C in continuous darkness. Six seedlings with nearly the same root length were then held in foam support floats that were suspended in 2.5-L plastic containers without Al3+ for acclimation to hydroponic conditions (0.5 mM CaCl2, pH 4.5, 16 h light/8 h dark). After 24 h of acclimation, the seedlings were photographed carefully using a camera (Nikon, COOLPIX A1000) to determine the main root lengths with a ruler beside them as scale. Then, the seedlings were transferred to solutions with or without AlCl3 (0.5 mM CaCl2, pH 4.5). The roots of the seedlings were photographed again after Al3+ exposure. To ensure the accuracy of the two kinds of measurements before and after Al3+ exposure, we marked the root at the initial position of the measurement. During the process of cultivation, the nutrient solution was aerated constantly with a flexible pipe connected with air pump.

The main root lengths were determined from the photographs using ImageJ software (National Institutes of Health, http://imagej.nih.gov/ij/). Root elongation was defined as the difference between the initial length before Al3+ treatment and the final length after Al3+ treatment. The root elongation under control (REC) and the root elongation under Al3+- stress (REA) were calculated, and the RRE was equal to REA/REC ×100% [37].

After Al3+ treatment, apical roots (0–2 cm) were excised by a scalpel, washed three times with 0.5 mM CaCl2 solution, and dried on filter paper. Then, six root tips for each line and the two parents were placed in a microcentrifuge tube (1.5 ml) containing 1.0 ml of 2 M HCl and extracted for 48 h with continuous shaking to release Al3+ from the soybean roots. The Al3+ levels in the extracts were determined by inductively coupled plasma-optical emission spectrometry (ICP-OES) (VARIAN 710-ES, America) [38].

Genotyping was carried out as previously described [35]. The soybean reference genome from Williams 82 [39] was used for read mapping be comparison with the tag sequence by SOAP software (The Beijing Genomics Institute, http://soap.genomics.org.cn/). Input data for SNP calling with realSFS was prepared by SAMtools [40]. RealSFS was used for SNP calling of every locus in the RIL population. The likelihoods of genotypes for each individual were integrated and extracted as candidate SNPs that were then filtered using the following criteria: 40 ≤ depth ≤ 2500, sites with a probability ≥ 95%. These highly reliable SNPs were used to obtain the genotypes of the parents and the RIL population. Moreover, the genotypes of all SNPs from the soybean genome were analyzed by the sliding window method and further used for each individual to generate bin information. Finally, a fine genetic map including 3,426 bin markers was constructed using MSTMap (http://alumni.cs.ucr.edu/yonghui/mstmap.html) and MapChart software (Wageningen University, https://www.wur.nl/en/show/Mapchart.htm) [41].

A high-density genetic map was constructed as previously described [35]. Composite interval mapping (CIM) was performed to detect QTLs using WinQTLCart software (North Carolina State University, http://statgen.ncsu.edu/qtlcart/WQTLCart.htm). The significant LOD threshold of 2.5 for QTLs was determined by a genome-wide permutation test with 1000 replications at the 5% level of significance. The analysis results also showed the effects of QTLs, the explanation rate of the phenotypic variation by QTLs and the interactions of QTLs. QTL mapping results were comprehensively compared to those published on Soybase (http://www.soybase.org/) [35].

The genes within all the QTL regions were listed by the Soybase website (http://www.soybase.org/). In addition, data from NCBI (https://www.ncbi.nlm.nih.gov/) and Phytozome (https://phytozome.jgi.doe.gov/pz/portal.html) were used to ascertain the conserved domains of the proteins and the possible functions of these domains. Specific primers for RT-PCR of these genes were designed using Primer Premier 5 software (PREMIER Biosoft, http://www.premierbiosoft.com/primerdesign/index.html).

The hydroponic conditions for cultivation of soybean seedlings were the same as those used for genotype analysis, as described in the section titles “Hydroponics and trait measurement”. Samples of the apical roots (0–2 cm) of the two parents were obtained and immediately frozen by using liquid nitrogen. Total RNA was extracted from the apical roots of seedlings grown under Al3+ stress or the control treatment using TRIzol reagent (TIANGEN, China). First-strand cDNA was synthesized using the PrimeScript™ RT Reagent Kit with gDNA Eraser (TAKARA, China) and used for further analysis of expression patterns for candidate genes.

The RT-PCR assay was carried out to analyze the expression of the genes in the apical roots from the two parent seedlings, with the soybean β-Tubulin gene as an internal reference, with the specific primers 5′-AACCTCCTCCTCATCGTACT-3′ and 5′-GACAGCATCAGCCATGTTCA-3′ [42]. The total volume of the PCR mixture was 20 μl, containing 1 μl of first-strand cDNA, 1 μl of each primer, 7 μl of ddH2O, and 10 μl of the mixture containing Taq DNA polymerase. The amplification reaction was performed as follows: predenaturation at 95°C for 3 min, followed by 30 cycles (for almost all genes; for β-Tubulin, 26 cycles were used) of 15 s at 95°C, 15 s at 54°C and 30 s at 72°C min, with a final extension for 5 min at 72°C. The PCR products were separated by agarose gel electrophoresis. Furthermore, qRT-PCR was further used to analyze the expression of the candidate genes. All PCRs were performed in 20-μl reactions consisting of 1 μl of cDNA, 0.8 μM each gene-specific primer and a mixture from the SYBR Green Supermix Kit (Takara, Japan). The reaction conditions were as follows: predenaturation at 94°C for 3 min, followed by 40 cycles of denaturation at 94°C for 10 s and renaturation at 54°C for 10 s; at the end of the reaction, the system was maintained at 95°C for 10 s, followed by lowering the temperature to 65°C for 5 s. The soybean Actin 3 gene [43] was used as an internal reference, with the forward primer 5′- GTGCACAATTGATGGACCAG-3′ and the reverse primer 5′-GCACCACCGGAGAGAAAATA-3′. Specific primers for RT-PCR and qRT-PCR of these genes were designed using Primer Premier 5 software (PREMIER Biosoft, http://www.premierbiosoft.com/primerdesign/index.html) (S1 and S2 Tables).

Analysis of variance (ANOVAs) was performed using SAS 9.4 by the general linear model (GLM) procedure with a logarithmic transform of data if necessary [44]. The broad-sense heritability (h2) of RRE and AAC was calculated according to Knapp et al [45]. Heritability was calculated using the formula shown as follows: h2 = σg2 / ((σe2 / n) + σg2), where σg2 denotes the genetic variance; σe2 denotes the error variance; and n denotes the replication number. The coefficient of variation was estimated as σg/μ, where μ represents the mean value. Phenotypic Pearson’s correlations were calculated using the ‘PROC CORR’ option of the SAS program between the two different traits [45]. Linear regression analysis was plotted using ‘MASS’ and ‘car’ package of R 3.5.4.

To explore the appropriate AlCl3 concentration and treatment time, two parents and five randomly selected lines were used to identify Al3+ tolerance characteristics. As shown in Fig 1, the RRE change trend was highly consistent among these lines. With increasing Al3+ concentrations, the RRE of each line decreased, showing strong inhibition of root elongation at high concentrations. Likewise, prolonged treatment of the lines led to decreased root elongation, so the 72-h treatment group exhibited the lowest RRE. In addition, the coefficient of variation (CV) was calculated to detect variations within each treatment. Comparative analysis showed that the condition with 25 μΜ [Al3+] and the 24-h treatment exhibited the highest CV (20.20%), which provided the greatest degree of dispersion among the five lines (Fig 1). Moreover, the two parents ZH 24 and HX 3 also showed the most significant difference under this condition (25 μΜ [Al3+], 24 h). Thus, 25 μΜ AlCl3 and 24 h of treatment were selected to obtain the widest separation in the RIL population.

https://doi.org/10.1371/journal.pone.0223674.g001

The ANOVA results demonstrated significant phenotypic differences among the RILs in RRE as well as AAC (P0.05), which was consistent with the characteristics of monogenic markers. There were 283 markers (account for 10.72%) showed separation distortion (P